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Terrestrial dark matter detection experiments probe the velocity-space distribution of dark matter 
particles in the vicinity of the Earth. We present a novel method, to be used in conjunction with 
standard cosmological simulations of hierarchical clustering, that allows one to extract a truly local 
velocity-space distribution in exquisite detail. Preliminary results suggest a new picture for this 
distribution which is decidedly non-Maxwellian but instead is characterized by randomly positioned 
peaks in velocity space. We discuss the implications of these results for both WIMP and axion 
detection experiments. 
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According to the cold dark matter (CDM) model of 
structure formation galaxies possess extended massive 
halos of non-baryonic particles. But while the CDM sce- 
nario can account for the rotation curves of spiral galax- 
ies and the morphology of structure on sub and super- 
galactic scales it sheds little light on the fundamental 
nature of dark matter. For this, we must turn to dark 
matter detection experiments. 

WIMPs (weakly interacting massive particles) and ax- 
ions are two dark matter candidates for which direct de- 
tection in a terrestrial apparatus is a distinct possibil- 
ity. Terrestrial detectors probe the distribution function 
(DF) of dark matter particles in the lab. The standard 
WIMP detection experiment attempts to measure the en- 
ergy deposited when a WIMP scatters off a nucleus in the 
detector Q. By contrast the direct detection scheme for 
axions relics on their coupling to photons: In a magnetic 
field, an axion can convert into a photon whose energy 
will be equal to the total (rest mass plus kinetic) energy 
of the axion 0]. 

The phase space DF, / (x, v) provides a complete de- 
scription of a collisionless system such as a dark matter 
halo. Formally / (x, v) d 3 xd 3 v is the number of parti- 
cles in the six-dimensional phase space element d 3 xd 3 v 
centered on the phase space point (x, v). A detector of 
volume V at position Xo probes the velocity space DF 
g Xo (v) = V -1 J v d 3 xf (x, v). Since detectors are very 
small, g Xo (v) = / (xo , v) to extremely high accuracy. 

The standard assumption for direct detection exper- 
iments is that g can be approximated by a Maxwellian 
with a cut-off at the escape speed of the Galaxy. Our aim 
is to explore the extent to which this assumption breaks 
down within the CDM scenario for structure formation. 
In the CDM model, dark halos are assembled through hi- 
erarchical clustering wherein small-scale objects collapse 
first and merge to yield systems of increasing size. Orig- 
inally, it was believed that the substructure one might 
have expected from hierarchical clustering was erased by 
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physical processes such as violent relaxation and tidal 
stripping However, recent high resolution simulations 
indicate that dark halos are in fact clumpy with a signifi- 
cant fraction of the mass bound in undigested lumps and 
tidal debris H0. 

The question is whether substructure in the density 
field of the halo translates into structure in the local 
velocity-space DF g. If so, there may be important im- 
plications for terrestrial detection experiments. (The im- 
plications of triaxiality, rotation, and a radial profile that 
differs from that of an isothermal sphere have been con- 
sidered by various authors 0, S H E3-) In the CDM 
scenario, the dark matter in the vicinity of the Earth is 
a superposition of bound clumps and tidal streams that 
were accreted by the Galaxy over its entire formation 
history. From the standpoint of direct detection experi- 
ments, particles from recent accretion events are the most 
interesting since they will have a bulk galactocentric ve- 
locity at the position of the Earth of order the escape 
speed of the Galaxy (~ 500 kms -1 ) or roughly twice the 
mean speed of the widely used Maxwellian model. More- 
over, even if most of the particles in these clumps have 
been stripped by the tidal field of the Galaxy, the tidal 
debris will have little time to phase mix and therefore 
be relatively coherent. Coherent streams of high veloc- 
ity particles may produce significant features in detection 
spectra even if they constitute a small fraction of the to- 
tal local dark matter density 0,0]- (See, however, [l3| 
where the authors reach a different conclusion.) 

The Maxwellian model may also be challenged by the 
following line of reasoning. The intrinsic velocity dis- 
persion of CDM is negligible and the distribution func- 
tion / (x, v) is most appropriately thought of as a three- 
dimensional sheet in six-dimensional phase space. Prior 
to the epoch of structure formation, this sheet is nearly 
perfectly flat: the velocities of particles at a particular 
point in space are given by the Hubble flow plus small 
perturbations. As structure formation proceeds, / re- 
tains its three-dimensional character but is "curled up" 
through the process of gravitational collapse, g (v) is 
the value of / on a three-dimensional surface that de- 
fines the detector, i.e., a surface that spans the three 
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velocity-space directions at the position of the detector. 
Since the intersection of two three-dimensional surfaces 
(one flat and the other highly contorted) is a collection of 
points, g will be a collection of discrete peaks in velocity 
space. These peaks will be randomly distributed though 
will likely exhibit clustering on different (velocity-space) 
scales. (For a general discussion of the dimensionality of 
phase-space structures as it relates to astronomical ob- 
servations, see Ref.[u|)- 

Fine structure of this type has been studied within the 
context of self-similar spherical infall models [l5j . These 
models assume that the initial density perturbation and 
angular momentum distribution are scale free and spher- 
ically symmetric. Given these assumptions, one can solve 
semi-analytically for the entire formation history of a 
halo and predict the positions of the velocity-space peaks 
which, in turn, may lead to a distinctive pattern of fea- 
tures in ener gy spectra obtained by axion and WIMP 
detectors |l0tll5j . 

Velocity-space structures in hierarchically formed ha- 
los can be studied using numerical simulations as follows: 
Consider a Milky Way-size simulated halo with N par- 
ticles and select a location that has the characteristics 
of the solar system (for example, a point 8.5 kpc from 
the center of mass of the halo) . Choose a volume V cen- 
tered on this point. The N' particles inside this volume 
provide a monte carlo realization of g(v). This proce- 
dure has been employed by various groups with some 
success However, even the largest simulations 

are unable to resolve velocity space structures deep in- 
side highly evolved halos. The limit on the velocity space 
resolution due to the finite number of particles in V is 
Av/cr ~ TV'- 1/3. Velocity space resolution is also limited 
by the finite size of V. Let L and $ be respectively the 
size of the halo and depth of its gravitational potential, 
L' be the size of V, and A<3> be the variation in <E> across 
V. The limit to the resolution in velocity space due to the 
finite size of V is Av/a ~ A$/$ - L'/L ~ (N'/N) 1/3 . 
We minimize Av/a by balancing the two effects (large 
V for particle number, small V for locality) to find the 
optimal choice N' ~ TV 1 / 2 . The implication is that for a 
simulation of N particles, only O (iV 1//2 ) are available to 
map the local velocity space distribution function. 

This paper introduces a new technique, to be used 
in conjunction with standard cosmological simulations, 
which allows one to map <?(v) at a single point within 
a dark halo. The method relies on test particles which, 
by design, reach the desired point (the position of the 
would-be detector) in the final timestep of the simulation. 
Through an iterative process, these test particles allow 
one to locate the points where the phase space sheet de- 
scribing the dark matter distribution intersects the phase 
space sheet describing the detector. Our method has a 
resolution in velocity space of Av/a ~ TV -1 / 3 as com- 
pared with Av/a ~ TV" 1 / 6 for the standard approach. 
Thus, with N = I0 6 , we can achieve a resolution that 
would have previously required 10 12 particles. 

The starting point for our algorithm is a simulation 



of the formation of a dark matter halo from standard 
CDM initial conditions. In the final frame of the simu- 
lation, we identify an appropriate location for a detec- 
tor and lay down a uniform grid (i.e., three-dimensional 
phase-space sheet) of massless test particles which spans 
velocity-space at that location. We next evolve both 
the test particles and original simulation particles back- 
ward in time to the initial frame of the simulation. The 
backward evolution is accomplished by simply reversing 
the timestep. During this "reverse run" the phase-space 
sheet defined by the test particles will fold and curl. At 
the initial time, we locate the intersection points of the 
test-particle phase-space sheet and the phase-space sheet 
that defines the initial dark matter DF. These intersec- 
tion points map directly into points in velocity space at 
the position of the detector (ie. the initial velocity-space 
distribution of the test particles). The value of the DF 
is calculated at each crossing point, and by Liouville's 
theorem, the DF of the test particles at the location of 
the detector can be determined. 

One may think of the test particles as providing the 
means to interpolate the DF in the region of the detec- 
tor. For the interpolation to make sense, the system 
must reverse properly, i.e., the positions of the simu- 
lation particles at the end of the backward integration 
should be within some small distance of their actual ini- 
tial positions. In essence, reversibility is equivalent to the 
condition that the dark matter DF, as described by the 
simulation particles, maintains its character as a three- 
dimensional structure in phase space. Chaos, which is 
almost certainly present in the orbits of the simulation 
particles, combined with round-off errors, has the poten- 
tial to spoil the reversibility of the time integration and 
destroy the integrity of the phase-space sheet. To apply 
our technique to a simulation, it must be verified that it 
does reverse properly. We have found that the simplest 
means of suppressing chaotic behavior and achieving re- 
versibility is to increase the softening length above what 
is normally used in cosmological simulations. 

At this point, the reader may suspect a swindle. If our 
method requires a larger softening length than is nor- 
mally used, how can we claim a substantial improvement 
over the standard N-body approach. Of course, chaos 
also affects the results of the standard method. We sus- 
pect that the ap prox imately Maxwellian DFs found in 
previous studies |l3t arise in part from the use of 
a finite volume to measure the velocity-space DF at a 
single point together with the fact that the dark matter 
phase-space sheet is not adequately modelled. Our in- 
terpolation scheme overcomes the finite volume problem 
and makes explicit the requirement that the phase space 
sheet be properly modelled. 

Typically we begin with the same number of test par- 
ticles as simulation particles so that the additional com- 
putation time is increased by a factor of less than 2. 
Higher resolution is achieved by an iterative procedure - 
regions surrounding intersection points are resampled by 
a finer grid of test particles and the reverse-run algorithm 
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is repeated. Our simulations use an implementation of 
Dehnen's O(N) algorithm modified to efficiently in- 
corporate test particles. Our technique however, is inde- 
pendent of the code used as long as it is time-reversible 
and allows test particles. The softening kernel is cho- 
sen such that the gravitational force goes to zero at zero 
separation but joins smoothly with Newtonian gravity at 
the softening length. In this manner, a test particle and 
a real massive particle with exactly the same phase-space 
coordinates will follow the same trajectory. 

As a demonstration of the power of our algorithm, we 
first examine the collapse of a spherically symmetric den- 
sity perturbation reminiscent of the self-similar collapse 
considered in [l5| and references therein. The pertur- 
bation is a monotonically decreasing function of radius 
so that the inner parts of the system collapse first. The 
simulation is run from t — to t = 12 (in units with 
Aftotai = 100 and G = 1) with approximately 33, 000 
particles and uses the same N-body code that will be 
used in subsequent, more realistic simulations. As the 
simulation is done in physical coordinates, all particles 
have an initial Hubble velocity and comoving softening 
length of 1.7 (Plummer equivalent) was used. 

The top panel of Figure^shows the final DF (radius r 
vs. radial velocity v r ). The fine structure (thin streams of 
particles) should be viewed as phase space sheets where 
we have projected out the 9 and <j> directions. In the ab- 
sence of spherical symmetry, the three-dimensional na- 
ture of the DF would still be present but not readily ap- 
parent in anr-t) r plot. Indeed, toward the center of the 
halo, where two-body scattering has destroyed spherical 
symmetry, the structures appear to dissolve. 

A detector at xq = (2,0,0) (horizontal line in the top 
panel of Figure ^) should measure 11 distinct velocity- 
space peaks. The middle panel of Figure 2] shows the 
results for the velocity space distribution function ob- 
tained using the standard N-body method. Two differ- 
ent volumes, one of radius Ar = 0.89 and containing 182 
particles (approximately N 1 / 2 ) and the other of radius 
Ar = 45 and containing 18 particles were used to gen- 
erate these results. With a large volume, one obtains 
a broad distribution and sees little evidence for discrete 
structures. A smaller volume picks out only a few of the 
peaks, and with very poor resolution. 

The bottom panel shows the results of the reverse-run 
method. All 11 peaks are detected. The plot is an actual 
histogram: the thickness of the peaks reflects the actual 
velocity-space resolution obtained by our method. Along 
the top of the plot are the positions of the peaks as mea- 
sured using the N-body approach but exploiting spherical 
symmetry, i.e., taking all particles close to r = 2. We see 
that the peaks from the reverse run method are precisely 
where they should be. 

The formation of a dark halo in the hierarchical clus- 
tering scenario is modelled by augmenting the spherically 
symmetric simulation considered above with a power-law 
spectrum (n = —2) of density and velocity perturbations 
(see, for example, reference |T^V The particle number 
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FIG. 1: Top: The r — v r diagram at the final time. Middle: 
The velocity distribution determined from the forward run 
with a large (solid) and small (dashed) volume at xo- Bottom: 
The reconstructed velocity distribution from our reverse-run 
technique. The dots along the top of the figure indicate the 
location of streams from the forward run when spherical sym- 
metry is applied. 



was increased to 5 x 10 5 and the softening set to 20 kpc 
Plummer equivalent. The simulation started at z = 44 
and evolved to the present in 1000 timesteps. The total 
system mass was 10 13 M Q . It should be noted that this 
simulation is only designed to qualitatively model the 
formation of the Milky Way so rigorous normalization 
of the power spectrum was not performed. Nonetheless, 
the resulting system consists of a massive central galaxy 
(M ~ 10 12 Mq) which partly formed through the merger 
of smaller objects and several dwarf satellites galaxies 
surrounding it. One can view the particle distribution 
and evolution of the simulation elsewhere [l^. Figure 
12 shows the speed distribution and a scatter plot in the 
two-dimensional phase space of speed and angle where 
the angle is measured with respect to some arbitrary di- 
rection in the halo. The latter may be of interest to 
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FIG. 2: Top panel: Speed distribution from the forward run 
with jV 1//2 particles (dotted) and from our reverse-run tech- 
nique (solid) . Lower panel: Distribution of angle between the 
local bulk motion and the crossing points. 



detectors with directional sensitivity such as the DRIFT 
detector |2(J. 

The top-left panel of Figure |3 gives an example of a 
recoil spectra that would be obtained by our fictitious 
detector. Since terrestrial detection experiments operate 
in the rest frame of the Earth we first shift the velocity 
distribution by 220 km s -1 . (We have not attempted to 
model the Galactic disk and therefore the direction of 
this shift is arbitrary.) In addition, we include the small 
velocity-space shift due to the Earth's motion around the 
Sun. It is this shift that gives rise to the famous seasonal 
modulation effect considered to be so important to ter- 
restrial detectors [2lJ. The seasonal modulation factor, 
S(E) = [Rate(June) - Rate(December)]/[Rate(June) + 
Rate(December)], (see Ref. 11]) is shown in the lower- 
left panel of the Figure. 

The recoil spectrum is characterized by a sequence of 
ste ps a s one goes up in energy. As noted in Refs. [ToL 
[Til Il2| peaks in the speed distribution translate to step 
functions in a recoil spectrum. 

Our ability to identify velocity space structure may be 
considerably better if dark matter is composed of axions. 
The standard axion detector measures the energy of the 
axion (kinetic plus rest mass) and therefore peaks in ve- 
locity space translate into peaks in the energy spectrum. 
In addition, since these experiments rely on the orienta- 
tion of the magnetic field, one will have directional sen- 
sitivity and therefore be able to map out the full <?(v). 
Figure |31 provides an example of an axion spectrum ob- 
tained by again shifting the velocity distribution to a 
terrestrial frame. 

The reverse- run algorithm provides the means to probe 
the local velocity-space structure of a simulated halo with 
a resolution far greater than is permitted through tradi- 



FIG. 3: Left Top: WIMP recoil energy spectra from the 
reverse-run (solid line) and a numerical Maxwellian distribu- 
tion designed to match the average properties of the reverse- 
run results (dotted line). Left Bottom: Seasonal modula- 
tion factor S(E). Right Top: Axion energy spectra from the 
reverse-run (solid line) and the same Maxwellian (dotted line). 
Right Bottom: Seasonal modulation factor S(E). 



tional methods. In short, our method allows us to obtain 
the velocity-space distribution at a given point without 
having to average over the surrounding region. 

Our results must be viewed as preliminary since our 
simulations focus a small region of the Universe and 
therefore do not take into account cosmological tidal 
fields. Moreover, no attempt is made to model disk for- 
mation. However, our method is completely general and 
may be applied to more realistic simulations of the Local 
Group or Local Supercluster. 

The preliminary results do suggest a new picture for 
the local velocity-space distribution function, one which 
is characterized by discrete peaks. This result is poten- 
tially of great importance to terrestrial detection exper- 
iments which rely on models of this distribution func- 
tion to develop search strategies and interpret experi- 
mental results. For WIMP searches, the deviations from 
Maxwellian and clumpy distributions are likely not sig- 
nificant for the current generation of detectors. However, 
current axion detectors have exceptional energy resolu- 
tion on the order of AE/(m a c 2 ) ~ 10" 11 22]. This res- 
olution is sufficient to be sensitive to seasonal and even 
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diurnal variations in the spectrum which can be used 
to fully resolve the velocity-space distribution of axions. 
Future generations of WIMP detectors with directional 
sensitivity, or accurate seasonal spectra may also yield 
detailed features of the dark matter distribution. 

We are in fact led to the tantalizing conclusion that 
if terrestrial experiments are able to detect dark matter, 



then the results will provide a wealth of information on 
the structure and formation history of the dark halo. 
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